#include <iostream>
#include "Types.hpp"
#include "Debug.hpp"
#include "ClWrapper.hpp"
#include "DeviceWrapper.hpp"
#include "FileReader.hpp"
#include "VTKPathWriter.hpp"
#include <algorithm>
#include <vector>
#include <fstream>
#include <string>
#include <cmath>
#include <iomanip>
#include <sys/time.h>
double timeStamp();


int main(int argc, char **argv)
{

	if(argc!=2){
		ERROR("Please Type the name of the input file");
		return 1;
	}
	std::string filename=argv[1];
	INFOP("");

	//TESTING READER
	FileReader reader=FileReader();

	//redistering keys
	reader.RegisterIntParameter(p_PART_OUT_FREQ,v_PART_OUT_FREQ);
	reader.RegisterIntParameter(p_VTK_OUT_FREQ,v_VTK_OUT_FREQ);
	reader.RegisterIntParameter(p_CL_WORKGROUP_SIZE,v_CL_WORKGROUP_SIZE);
	reader.RegisterIntParameter(p_CLWORKGROUP_1DSIZE,v_CLWORKGROUP_1DSIZE);
	reader.RegisterIntParameter(p_CLWORKGROUP_3DSIZE_X,v_CLWORKGROUP_3DSIZE_X);
	reader.RegisterIntParameter(p_CLWORKGROUP_3DSIZE_Y,v_CLWORKGROUP_3DSIZE_Y);
	reader.RegisterIntParameter(p_CLWORKGROUP_3DSIZE_Z,v_CLWORKGROUP_3DSIZE_Z);
	reader.RegisterIntParameter(p_XN,v_XN);
	reader.RegisterIntParameter(p_YN,v_YN);
	reader.RegisterIntParameter(p_ZN,v_ZN);
	reader.RegisterIntParameter(p_REFLECTX,v_REFLECTX);
	reader.RegisterIntParameter(p_REFLECTY,v_REFLECTY);
	reader.RegisterIntParameter(p_REFLECTZ,v_REFLECTZ);

	reader.RegisterDoubleParameter(p_TIME_STEP,v_TIME_STEP);
	reader.RegisterDoubleParameter(p_END_TIME,v_END_TIME);
	reader.RegisterDoubleParameter(p_EPSILON,v_EPSILON);
	reader.RegisterDoubleParameter(p_SIGMA,v_SIGMA);
	reader.RegisterDoubleParameter(p_X_MIN,v_X_MIN);
	reader.RegisterDoubleParameter(p_X_MAX,v_X_MAX);
	reader.RegisterDoubleParameter(p_Y_MIN,v_Y_MIN);
	reader.RegisterDoubleParameter(p_Y_MAX,v_Y_MAX);
	reader.RegisterDoubleParameter(p_Z_MIN,v_Z_MIN);
	reader.RegisterDoubleParameter(p_Z_MAX,v_Z_MAX);
	reader.RegisterDoubleParameter(p_R_CUT,v_R_CUT);
	reader.RegisterDoubleParameter(p_KS,v_KS);
	reader.RegisterDoubleParameter(p_KD,v_KD);
	reader.RegisterDoubleParameter(p_GX,v_GX);
	reader.RegisterDoubleParameter(p_GY,v_GY);
	reader.RegisterDoubleParameter(p_GZ,v_GZ);

	reader.RegisterStringParameter(p_INPUT_FILE,v_INPUT_FILE);
	reader.RegisterStringParameter(p_PART_OUT_NAME,v_PART_OUT_NAME);
	reader.RegisterStringParameter(p_VTK_OUT_NAME,v_VTK_OUT_NAME);

	//Parsing parameters
	reader.ReadFile(filename);
	//cl_real epsilon=reader.GetDoubleParameter(p_EPSILON);
	//cl_real sigma=reader.GetDoubleParameter(p_SIGMA);
	//cl_real sigma6=pow(sigma,6);
	cl_uint WORK_GROUP_SIZE=reader.GetIntParameter(p_CL_WORKGROUP_SIZE);
	cl_real delta_t=reader.GetDoubleParameter(p_TIME_STEP);
	cl_real end_t=reader.GetDoubleParameter(p_END_TIME);
	int vtkFreq = reader.GetIntParameter(p_VTK_OUT_FREQ);
	int partfreq = reader.GetIntParameter(p_PART_OUT_FREQ);
	int reflectx = reader.GetIntParameter(p_REFLECTX);
	int reflecty = reader.GetIntParameter(p_REFLECTY);
	int reflectz = reader.GetIntParameter(p_REFLECTZ);
	cl_uint nx =  reader.GetIntParameter(p_XN);
	cl_uint ny =  reader.GetIntParameter(p_YN);
	cl_uint nz =  reader.GetIntParameter(p_ZN);
	cl_real x_min =reader.GetDoubleParameter(p_X_MIN);
	cl_real x_max =reader.GetDoubleParameter(p_X_MAX);
	cl_real y_min =reader.GetDoubleParameter(p_Y_MIN);
	cl_real y_max =reader.GetDoubleParameter(p_Y_MAX);
	cl_real z_min =reader.GetDoubleParameter(p_Z_MIN);
	cl_real z_max =reader.GetDoubleParameter(p_Z_MAX);
	cl_real len_x= (x_max-x_min)/nx;
	cl_real len_y=  (y_max-y_min)/ny;
	cl_real len_z=  (z_max-z_min)/nz;
	//cl_real cut_r= reader.GetDoubleParameter(p_R_CUT);
	cl_real ks = reader.GetDoubleParameter(p_KS);
	cl_real kd = reader.GetDoubleParameter(p_KD);
	cl_real min1=std::min(len_x,len_y);
	cl_real min_l=std::min(len_z,min1);

	cl_real Gx = reader.GetDoubleParameter(p_GX);
	cl_real Gy = reader.GetDoubleParameter(p_GY);
	cl_real Gz = reader.GetDoubleParameter(p_GZ);


	// Asserting the Cut Off Radius to be Smaller than the Minimum of the Cell Sizes
	//if(!(cut_r<=len_x && cut_r<=len_y && cut_r<=len_z)){
	//std::cerr<<"Cut-Off Radius is greater than the Cell Size- Exiting";
	//return EXIT_FAILURE;
	//}

	std::string part_out_name= reader.GetStringParameter(p_PART_OUT_NAME);
	std::string vtk_out_name= reader.GetStringParameter(p_VTK_OUT_NAME);

	// Getting information about the Platforms
	ClWrapper platofrmCL;
	std::cout << "Available Platforms" << std::endl;
	std::cout << platofrmCL;

	// Get the Platform from the USER
	cl::Platform pl = platofrmCL.getPlatform();	
	std::cout << "Using the Following Platform: " << std::endl;
	ClWrapper platform_sel(pl);
	std::cout << platform_sel;

	// Obtain information about the Devices present in the Selected Platform
	std::cout << "Available Devices" << std::endl;
	DeviceWrapper deviceCL(pl);
	std::cout << deviceCL;

	cl::Device device;      
	cl::Context context; 
	cl::CommandQueue queue; 
	cl::Program::Sources sources;
	cl::Program program;

	//Set the Device, Context, Program, Commandqueue, sources. Can be modified later to make number of Programs instead of one program.                                                                  
	deviceCL.getDeviceProgram("a2_kernels.cl", device, context, queue, sources, program);
	std::cout << "Using the Following Device: " << std::endl;
	DeviceWrapper deviceSel(device);
	std::cout << deviceSel;


	//reader.PrintParameters(std::cout);
	cl_uint N;
	cl_uint NCells=nx*ny*nz;
	cl_real* masses;
	cl_real* radius;
	cl_real4* positions;
	cl_real4* velocities;
	cl_int* cells=new cl_int[NCells];

	std::string fileParti=reader.GetStringParameter(p_INPUT_FILE);
	reader.readParticleFile(N,masses,radius,positions,velocities,fileParti);
	cl_int* particle_ids=new cl_int[N]; // Was defined before reading N from the file- Corrected
	cl_real4* oldForces=new cl_real4[N];
	cl_real4* forces=new cl_real4[N];
	cl_uint GLOBAL_CELLS_SIZE=ceil(NCells/(float)WORK_GROUP_SIZE)*WORK_GROUP_SIZE;
	cl_uint GLOBAL_PARTICLES_SIZE=ceil(N/(float)WORK_GROUP_SIZE)*WORK_GROUP_SIZE;

	for(size_t i=0;i<N;++i)
	{
		if(radius[i]>min_l*0.5)
		{
			std::cout<<"ERROR: Particle "<<i<<" has a radius " << radius[i] << "> L/2\n";
			return 1;
		}
	}

	//reader.PrintParameters(std::cout);
	//Print the values from the file
	/*for(cl_uint i=0;i<N;++i)
	{
	std::cout<<masses[i]<<"  "<<radius[i]<<"  "<<positions[i].s[0]<<"  "<<positions[i].s[1]<<"  "<<positions[i].s[2]<<"  "
	<<velocities[i].s[0]<<"  "<<velocities[i].s[1]<<"  "<<velocities[i].s[2]<<"\n";
	}*/


	try
	{


		//Positions buffer
		cl::Buffer posBuffer=cl::Buffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR ,sizeof(cl_real4) * N, &positions[0]);
		//Velocities buffer
		cl::Buffer velBuffer=cl::Buffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR ,sizeof(cl_real4) * N, &velocities[0]);
		//radius buffer
		cl::Buffer radiBuffer=cl::Buffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR ,sizeof(cl_real) * N, &radius[0]);

		cl::Buffer cellsBuffer=cl::Buffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR ,sizeof(cl_int) * NCells, &cells[0]);

		cl::Buffer part_id_Buffer=cl::Buffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR ,sizeof(cl_int) * N, &particle_ids[0]);

		//(1/2m) buffer of masses
		cl::Buffer inv2MassBuffer = cl::Buffer(context, CL_MEM_ALLOC_HOST_PTR, sizeof(cl_real) * N);
		cl_real* inv_2masses = (cl_real*)queue.enqueueMapBuffer(inv2MassBuffer, true, CL_MAP_WRITE, 0,N * sizeof(cl_real));
		for(cl_uint i=0;i<N;++i)
			inv_2masses[i]=1/(2*masses[i]);
		queue.enqueueUnmapMemObject(inv2MassBuffer, (void*)inv_2masses);

		//Force buffers
		cl::Buffer forceBuffer = cl::Buffer(context, CL_MEM_ALLOC_HOST_PTR, sizeof(cl_real4) * N);
		cl::Buffer oldForceBuffer = cl::Buffer(context, CL_MEM_ALLOC_HOST_PTR, sizeof(cl_real4) * N);


		// masses*G buffer of Gravitational Force
		cl::Buffer gravForceBuffer = cl::Buffer(context, CL_MEM_ALLOC_HOST_PTR, sizeof(cl_real4) * N);
		cl_real4* gravForce = (cl_real4*)queue.enqueueMapBuffer(gravForceBuffer, true, CL_MAP_WRITE, 0,N * sizeof(cl_real4));

		for(cl_uint i=0;i<N;++i)
		{
			if (masses[i]!=std::numeric_limits<cl_real>::infinity())
			{
				gravForce[i].s[0]=masses[i]*Gx;
				gravForce[i].s[1]=masses[i]*Gy;
				gravForce[i].s[2]=masses[i]*Gz;
				}
			else
			{
				gravForce[i].s[0]=0.0;
				gravForce[i].s[1]=0.0;
				gravForce[i].s[2]=0.0;

			}
			gravForce[i].s[3]=0.0;			
		}
		queue.enqueueUnmapMemObject(gravForceBuffer, (void*)gravForce);


		//compute the forces
		cl::Kernel kCalcForces = cl::Kernel(program, "forceCalculation");
		kCalcForces.setArg(0, posBuffer);
		kCalcForces.setArg(1, forceBuffer);
		kCalcForces.setArg(2, cl::__local(sizeof(cl_real4) * WORK_GROUP_SIZE));
		kCalcForces.setArg(3, cellsBuffer);
		kCalcForces.setArg(4, part_id_Buffer);
		kCalcForces.setArg(5, radiBuffer);
		kCalcForces.setArg(6, velBuffer);
		kCalcForces.setArg(7, N);
		kCalcForces.setArg(8, nx);
		kCalcForces.setArg(9, ny);
		kCalcForces.setArg(10, nz);
		kCalcForces.setArg(11, len_x);
		kCalcForces.setArg(12, len_y);
		kCalcForces.setArg(13, len_z);
		kCalcForces.setArg(14, ks);
		kCalcForces.setArg(15, kd);
		kCalcForces.setArg(16, reflectx);
		kCalcForces.setArg(17, reflecty);
		kCalcForces.setArg(18, reflectz);
		kCalcForces.setArg(19, x_min);
		kCalcForces.setArg(20, y_min);
		kCalcForces.setArg(21, z_min);
		kCalcForces.setArg(22, x_max);
		kCalcForces.setArg(23, y_max);
		kCalcForces.setArg(24, z_max);
		kCalcForces.setArg(25, gravForceBuffer);

		//Update positions and copy forces to oldForces
		cl::Kernel kUpdatePositions = cl::Kernel(program, "updatePositions");
		kUpdatePositions.setArg(0, posBuffer);
		kUpdatePositions.setArg(1, oldForceBuffer);
		kUpdatePositions.setArg(2, forceBuffer);
		kUpdatePositions.setArg(3, velBuffer);
		kUpdatePositions.setArg(4, inv2MassBuffer);
		kUpdatePositions.setArg(5, delta_t);
		kUpdatePositions.setArg(6, N);

		//Update velocities
		cl::Kernel kUpdateVelocities = cl::Kernel(program, "updateVelocities");
		kUpdateVelocities.setArg(0, velBuffer);
		kUpdateVelocities.setArg(1, oldForceBuffer);
		kUpdateVelocities.setArg(2, forceBuffer);
		kUpdateVelocities.setArg(3, inv2MassBuffer);
		kUpdateVelocities.setArg(4, delta_t);
		kUpdateVelocities.setArg(5, N);

		//set Cells to default empty value
		cl::Kernel kSetCellsToMinus = cl::Kernel(program, "setCellsToMinus");
		kSetCellsToMinus.setArg(0,cellsBuffer);
		kSetCellsToMinus.setArg(1,NCells);

		//Assemble Cells
		cl::Kernel kAssembleCells = cl::Kernel(program, "assembleCells");
		kAssembleCells.setArg(0,posBuffer);
		kAssembleCells.setArg(1,part_id_Buffer);
		kAssembleCells.setArg(2,cellsBuffer);
		kAssembleCells.setArg(3,N);
		kAssembleCells.setArg(4,nx);
		kAssembleCells.setArg(5,ny);
		kAssembleCells.setArg(6,nz);
		kAssembleCells.setArg(7,len_x);
		kAssembleCells.setArg(8,len_y);
		kAssembleCells.setArg(9,len_z);

		cl::Kernel kSetPeriodicity = cl::Kernel(program, "setPeriodicity");
		kSetPeriodicity.setArg(0,posBuffer);
		kSetPeriodicity.setArg(1,N);
		kSetPeriodicity.setArg(2,x_min);
		kSetPeriodicity.setArg(3,x_max);
		kSetPeriodicity.setArg(4,y_min);
		kSetPeriodicity.setArg(5,y_max);
		kSetPeriodicity.setArg(6,z_min);
		kSetPeriodicity.setArg(7,z_max);
		kSetPeriodicity.setArg(8, reflectx);
		kSetPeriodicity.setArg(9, reflecty);
		kSetPeriodicity.setArg(10, reflectz);


		/**************SIMULATION******************/

		//Initialize the particles ids and the cells to -1
		queue.enqueueNDRangeKernel(kSetCellsToMinus, cl::NullRange, cl::NDRange(GLOBAL_CELLS_SIZE), cl::NDRange(WORK_GROUP_SIZE));
		//	queue.enqueueNDRangeKernel(kSortParticleIDs,cl::NullRange, cl::NDRange(GLOBAL_PARTICLES_SIZE), cl::NDRange(WORK_GROUP_SIZE));
		queue.enqueueReadBuffer(cellsBuffer, CL_TRUE, 0, sizeof(cl_int) * NCells, cells);
		queue.enqueueReadBuffer(part_id_Buffer, CL_TRUE, 0, sizeof(cl_int) * N, particle_ids);

		queue.enqueueNDRangeKernel(kAssembleCells, cl::NullRange, cl::NDRange(GLOBAL_PARTICLES_SIZE), cl::NDRange(WORK_GROUP_SIZE));
		queue.enqueueReadBuffer(cellsBuffer, CL_TRUE, 0, sizeof(cl_int) * NCells, cells);
		queue.enqueueReadBuffer(part_id_Buffer, CL_TRUE, 0, sizeof(cl_int) * N, particle_ids);

		int counterp=0;
		for(cl_uint idC=0;idC<NCells;++idC){
			int id=cells[idC];
			if(id!=-1){
				std::cout<<"CELL "<<idC<<" Particles: ";
				while(id!=-1)
				{
					++counterp;
					std::cout<< id <<"  ";
					id = particle_ids[id];
				}
				std::cout<<std::endl;
			}
		}
		//Compute the initial forces
		queue.enqueueNDRangeKernel(kCalcForces, cl::NullRange, cl::NDRange(GLOBAL_PARTICLES_SIZE), cl::NDRange(WORK_GROUP_SIZE));
		int counter=0,n_print=0,p_print=0;
		typedef std::numeric_limits< double > dbl;
		
		for(cl_real t=0;t<end_t;t+=delta_t)
		{
			if(counter%partfreq==0)			
			{
				queue.enqueueReadBuffer(posBuffer, CL_TRUE, 0, sizeof(cl_real4) * N, positions);
				queue.enqueueReadBuffer(velBuffer, CL_TRUE, 0, sizeof(cl_real4) * N, velocities);
				//std::string filename="vtkFiles//";
				std::string filename=part_out_name;
				filename+=std::to_string(p_print);
				filename += ".out";
				std::ofstream mStream;
				mStream.open(filename);
				mStream << N<<"\n";

				for(cl_uint i(0);i<N;++i){
					mStream<<std::fixed;
					mStream<< std::setprecision(6) << masses[i] <<" "<< radius[i]<<" "<< positions[i].s[0]<<" "<<positions[i].s[1]<<" "<<positions[i].s[2]<<" ";					
					mStream<<velocities[i].s[0]<<" "<<velocities[i].s[1]<<" "<<velocities[i].s[2]<<"\n";
				}
				++p_print;
				mStream.close();			  
			}

			if(counter%vtkFreq==0)
			{
				queue.enqueueReadBuffer(posBuffer, CL_TRUE, 0, sizeof(cl_real4) * N, positions);
				queue.enqueueReadBuffer(velBuffer, CL_TRUE, 0, sizeof(cl_real4) * N, velocities);
				//std::string filename1="vtkFiles//";
				std::string filename1=vtk_out_name;
				filename1+=std::to_string(n_print);
				VTKPathWriter writer(filename1,N);
				writer.writePoints(positions);
				writer.writeScalars(masses,"m");
				writer.writeScalars(radius,"r");
				writer.writeVector(velocities,"v");
				++n_print;
			}
			++counter;
			std::cout<<"time\t" << t<<std::endl;
			//	std::cout<<"Update_POSITIONS"<<std::endl;
			queue.enqueueNDRangeKernel(kUpdatePositions, cl::NullRange, cl::NDRange(GLOBAL_PARTICLES_SIZE), cl::NDRange(WORK_GROUP_SIZE));
			//std::cout<<"SET_PERIODICITY\n";
			queue.enqueueNDRangeKernel(kSetPeriodicity, cl::NullRange, cl::NDRange(GLOBAL_PARTICLES_SIZE), cl::NDRange(WORK_GROUP_SIZE));
			//std::cout<<"SET_MINUS"<<std::endl;
			queue.enqueueNDRangeKernel(kSetCellsToMinus, cl::NullRange, cl::NDRange(GLOBAL_CELLS_SIZE), cl::NDRange(WORK_GROUP_SIZE));
			//	std::cout<<"SORT_PARTICLES"<<std::endl;
			//	queue.enqueueNDRangeKernel(kSortParticleIDs,cl::NullRange, cl::NDRange(GLOBAL_PARTICLES_SIZE), cl::NDRange(WORK_GROUP_SIZE));
			//	std::cout<<"ASSEMBLE_CELLS"<<std::endl;
			queue.enqueueNDRangeKernel(kAssembleCells, cl::NullRange, cl::NDRange(GLOBAL_PARTICLES_SIZE), cl::NDRange(WORK_GROUP_SIZE));
			//	std::cout<<"CALC_FORCES"<<std::endl;
			queue.enqueueNDRangeKernel(kCalcForces, cl::NullRange, cl::NDRange(GLOBAL_PARTICLES_SIZE), cl::NDRange(WORK_GROUP_SIZE));
			//	std::cout<<"UPDATE_VELOCITIES"<<std::endl;
			queue.enqueueNDRangeKernel(kUpdateVelocities, cl::NullRange, cl::NDRange(GLOBAL_PARTICLES_SIZE), cl::NDRange(WORK_GROUP_SIZE));
			//	std::cout<<"NEXT_\n";
		}
		queue.finish();
	}
	catch(cl::Error er) 
	{
		printf("ERROR: %s(%d)\n", er.what(), er.err());
	}

	delete [] masses;
	delete [] positions;
	delete [] velocities;
	delete [] forces;
	delete [] oldForces;
	delete [] cells;
	delete [] particle_ids;

	return 0;
}

double timeStamp()
{
	struct timeval tp;
	gettimeofday(&tp, NULL);
	return((double)(tp.tv_sec + tp.tv_usec/1000000.0));
}

